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Abstract 

Chemical reactions on dust grains are of crucial importance in interstellar chemistry because they 
produce molecular hydrogen and various organic molecules. Due to the submicron size of the grains 
and the low flux, the surface populations of reactive species are small and strongly fluctuate. Under 
these conditions rate equations fail and the master equation is needed for modeling these reactions. 
However, the number of equations in the master equation grows exponentially with the number 
of reactive species, severely limiting its feasibility. Here we present a method which dramatically 
reduces the number of equations, thus enabling the incorporation of the master equation in models 
of interstellar chemistry. 

PACS numbers: 05.10.-a,82.65.+r,98.58.-w 
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The chemistry of interstellar clouds consists of reactions taking place in the gas phase as 
well as on the surfaces of dust grains 11. Surface reactions include the formation of molecular 

nn 

hydrogen |2J, |3( and reaction networks producing various organic molecules. The simulations 
of grain-surface chemistry are typically done using rate equation models J, 0, 0, B, Isj] . These 
models consist of coupled ordinary differential equations that provide the time derivatives of 
the populations of the reactive species on the grains. Rate equations are highly efficient for 
the simulation of reactions on macroscopic surfaces. However, in the limit of small grains 



under low flux, rate equations fail because they ignore the discrete nature of the populations 



of reactive species and their fluctuations 
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Recently, a master equation approach for the simulation of reaction networks on small 



grains was proposed \14L Il5l |. It takes into account the discreteness and fluctuations in 
the surface populations of reactive species and provides accurate results for the reaction 
rates. For example, in the case of hydrogen recombination, its dynamical variables are 
the probabilities P{N) that there are N hydrogen atoms on a grain. The time derivatives 
P(N), N = 0,1,2,..., are expressed in terms of the adsorption, reaction and desorption 
terms. The master equation was applied to the study of reaction networks involving multiple 



16|, Il2|- The master equation can 



species that appear in models of interstellar chemistry 

be solved either by direct integration jl^l is| or by using a Monte Carlo (MC) method [18]. 
A significant advantage of direct integration over the MC approach is that the equations 
can be easily coupled to the rate equations of gas-phase chemistery. However, the number 
of coupled equations increases exponentially with the number of reactive species, making 
direct integration infeasible for complex reaction networks of multiple species |lfi. 17]. 

In this paper we introduce the multi-plane method, in which the number of equations 
is dramatically reduced, thus enabling the incorporation of the master equation in models 
of interstellar chemistry. The multi-plane method is tested by comparing its results to 
those obtained from the complete master equation set, showing excellent agreement. To 
demonstrate the method we consider a reaction network that involves three reactive species: 
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For simplicity we denote the reactive 



H and O atoms and OH molecules 
species by Xi = H, X 2 = O, X 3 = OH, and the resulting non- reactive species by X 4 = H 2 , 
X5 = O2, Xq = H 2 0. The reactions that take place in this network include H + O — > OH 
{X x + X 2 -> X 3 ), H + H -> H 2 {X x + X x -> X 4 ), O + O -> 2 (X 2 + X 2 -> X 5 ), and H + 
OH -> H 2 (X 1 + X 3 — X 6 ). 



2 



Consider a small spherical grain of diameter d, exposed to fluxes of H and atoms 
and OH molecules. The cross-section of the grain is a = ird 2 /4 and its surface area is 
nd 2 . The density of adsorption sites on the surface is denoted by s (sites cm -2 ). Thus, 
the number of adsorption sites on the grain is S = 7rd 2 s. The desorption rates of atomic 
and molecular species on the grain are given by Wi = v ■ exp[— P^^/Zc^T], where v is the 
attempt rate (standardly taken to be 10 12 s _1 ), E\[i) is the activation energy barrier for 
desorption of specie Xi and T (K) is the surface temperature. The hopping rate of adsorbed 
atoms between adjacent sites on the surface is a% = v ■ exp[— E (i)/kBT], where E (i) is 
the activation energy barrier for hopping of Xi atoms (or molecules). Here we assume that 
diffusion occurs only by thermal hopping, in agreement with experimental results For 
small grains, it is convenient to replace the hopping rate <2j (hops s -1 ) by the sweeping rate 
M = cti/S, which is approximately the inverse of the time it takes for an Xi atom to visit 
nearly all the adsorption sites on the grain surface (up to a logarithmic correction [2^). 

In the reaction network described above, the dynamical variables of the master equation 
are the probabilities P(Ni, N 2 , N 3 ) of having a population that includes Ni atoms of specie 
Xi on the grain. Only the reactive species are included in the master equation, from which 
the production rates of the non-reactive species can be obtained. The master equation for 
the H, O and OH system takes the form 

P(N U N 2 , N 3 ) = ]T P [P(., Ni-1,..)- P(N U N 2 , N 3 )] 

i=l 

+ m + 1)P(., Ni + l,..)- N t P(N u N 2 , N 3 )] 

i=i 
2 

+ A [{Ni + 2)(Ni + 1)P(.., Ni + 2, ..) - Ni(Ni - 1)P(N U N 2 , N 3 )) 
i=i 

+(A + A 2 ) [(JVi + 1)(N 2 + l)P(iV 1 + 1, N 2 + 1, N 3 - 1) - mNzPim, N 2 , N 3 )} 
+(A + A 3 ) [(JVi + 1)(N 3 + l)P(iV 1 + 1, N 2 , N 3 + l)- N 1 N 3 P(Ni, N 2 , N 3 )} . (1) 

The terms in the first sum describe the incoming flux, where Pj (atoms s -1 ) is the flux per 
grain of the specie Xi. The probability P(.., N, ..) increases due to adsorption of an Xi 
atom on grains that already have Ni — 1 such atoms on their surfaces, and decreases due 
to adsorption on grains that have Ni atoms. Similarly, the second sum describes the effect 
of desorption. The third sum describes the effect of diffusion mediated reactions between 
two atoms of the same specie and the last two terms account for reactions between different 
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species. Each reaction rate is proportional to the number of pairs of atoms of the two species 
involved, and to the sum of their sweeping rates. The average population size of the X{ specie 
on the grain is (Ni) = En,,n 2 ,n 3 N { P(N x , N 2 , N 3 ), where N { = 0, 1, 2 . . ., and % = 1, 2 or 
3. The production rate per grain R(X k ) (molecules s _1 ) of X k molecules produced by the 
reaction Xi+X d -> X k is given by R(X k ) = {A i + A j )(N i Nj), or by R(X k ) = ^(iV^AW)) 
in case that i = j. 

In numerical simulations the master equation must be truncated in order to keep the 
number of equations finite. A convenient way to achieve this is to assign upper cutoff's 
A^ max , i = 1, . . . , J on the population sizes, where J is the number of reactive species. The 
number of coupled equations is thus 

N E = f[(N™ ax + l). (2) 

i=l 

The truncated master equation is valid if the probability to have larger populations beyond 
the cutoffs is vanishingly small. However, the number of equations, Ne, grows exponentially 
as the number of reactive species increases. This severely limits the applicability of the 
master equation to interstellar chemistry. To exemplify the magnitude of this difficulty, 
note that the upper cutoffs must satisfy iVf iax > 1 for the J' species that do not react 
with themselves and j\T max > 2 for the J" species that do react with themselves. Therefore, 
N E > 2 J '3 J " 0- 

The chemistry taking place on interstellar dust grains is dominated by hydrogen atoms 
for three reasons. First, hydrogen is the most abundant specie; Second, it hops on the 
surface much faster than other species; Third, it is the most reactive specie, namely it reacts 
with many other species which do not react with each other. If two species such as X 2 and 
X 3 do not react with each other it is not crucial to maintain the correlation between their 
population sizes on a grain. Therefore, the probability distribution of the population sizes 
can be approximated by 

P(JV 1? N 2 , N 3 ) = P(N 1 )P(N 2 /N 1 )P(N 3 /N l ), (3) 

where P(Ni/Ni) is the conditional probability for a population size Ni of specie X it given 
that the population size of Xi atoms is N%. 

The multi-plane method is derived as follows: inserting Eq. © into the master equation 
(HJ), we trace over N 3 , using the fact that Ev 3 P(iV 3 /iVi) = 1 and Ejv 3 P(iV 3 /iVi) = 0. We 

4 



obtain the following set of equations: 

P(Aq, A 2 ) = ]T P [P(.., Ni-1,..)- P(N U N 2 )] 
i=i 

2 

+ E W i i( N i + Ni + l,..)- NiP(N u A 2 )] 

i=i 

2 

+ Y,Ai [(Ni + 2)(Ni + 1)P(.., N t + 2, ..) - A,(A, - 1)P(A 1; A 2 , N 2 )} 
i=i 

+{A X + A 2 ) [{Nt + 1)(N 2 + l)P(Ni + 1, N 2 + 1) - N 1 N 2 P(N 1 , N 2 )} 

+(Ax + A 3 ) [(JVj + l)P(JVj + 1, Nz)(N 3 ) Nl+1 - NtPi^, N 2 )(N 3 ) Nl ] , (4) 

where (N 3 ) Nl = J2n s N 3 P(N 3 /Ni). A similar procedure for tracing over A 2 yields 

P{N 1} N 3 ) = £ P [ p (-, ^ - 1, ..) - P(N U N 3 )} 

i=l,3 

+ £ [(^ + ^ + 1, ••) - ^p(a 15 ^ 3 )] 

i=l,3 

+A 1 [(N, + 2)(Aq + l)P(Ai + 2, A 3 ) - Ai(A r i - 1)P(A 1; A 3 )] 

+ (A + A 2 ) [(A x + 1)P(A X + 1, A 3 - 1) (i\T 2 ) jvi+i - A 1 P(A 1 , A 3 )(A 2 ) Wl ] 

+ A 3 ) [(A"i + 1)(A 3 + 1)P(A X + 1, A 3 + 1) - A 1 A 3 P(A 1 , A 3 )] . (5) 

The production rates of the non-reactive species X 4 , X 5 and A 6 are given by 
R(X t+3 ) = Ai J2n!,n 2 Ni(Ni - 1)P(A 1; A 2 ), where % = 1,2, and P(A 6 ) = (A x + 
A 3 ) J2ni n 3 NiNaP(Ni, N 3 ). The desorption rate of the reactive specie A 3 is given by 
R(X 3 ) = W 3 J2n!,n 3 N 3 P(Ni, A 3 ). In simulations using the multi-plane method, one can 
choose any initial condition that can be expressed by Eq. @. A convenient choice is of an 
empty grain, namely, P(A X = 0, A 2 = 0) = 1 in Eq. (gj) and P(A X = 0, A 3 = 0) = 1 in Eq. 
(|SJ), where all other probabilities vanish. The multi-plane method requires setting cutoffs, 
A™ ax , i — 1, 2, 3, where the same value of N™ ax is used in Eqs. (@J and (jSJ). 

The simulations presented here were done for spherical grains on which the density of 
adsorption sites is s = 5 x 10 13 (sites cm -2 ), at a grain temperature of T = 10 K. The 
activation energies for diffusion and desorption of the reactive species H, O and OH were 
taken as E (l) = 22, E X {1) = 32, E (2) = 25, E x {2) = 32 and E (3) = 28, Pi(3) = 35 meV, 
respectively. The parameters for hy drog en are rounded values in the range of experimental 
results on silicate and ice surfaces [l9|. For the other species no concrete experimental 
results are available and the chosen values reflect the tendency of heavier species to be more 



5 



strongly bound to the surface. The flux of H atoms on a grain of S sites was taken as 
F\ = 2.75 x lO^S 1 (s _1 ), and the flux of O atoms was F 2 = O.OlFi. These parameters are 
suitable for dense molecular clouds where such reaction networks are likely to take place. For 
simplicity, the OH flux was neglected, taking F 3 = 0. In Fig. 1 we present the production 
rates of H 2 (circles), 2 (squares), H 2 (triangles) and the desorption rate of OH (x) per 
grain vs. grain size. The results of the multi-plane method (symbols) and the complete 
master equation (solid lines), both obtained by direct numerical integration, are in excellent 
agreement. The results of the rate equations, that for small grains deviate significantly from 
the master equation results are also shown (dashed line). 

In a complex reaction network of J reactive species, Xj, i = 1, . . . , J, the distribution 
of surface populations is given by P(Ni, . . . , Nj). If the only reactions are of the form 
X\ + Xi — > X p and Xi + Xj — > X q (no matter if X p and X q are reactive or not), then the 
system can be reduced to J — 1 sets of equations for P(Ni, Nj), j = 2, . . . , J. Each of these 
sets is obtained by tracing over all the populations N except for Ni and Nj. The number 
of equations is reduced from the exponential form of Eq. (J2J) to 

j 

N E = (iVT* + 1) £(W? MX + 1), (6) 

namely to a linear dependence on the number of reactive species. In case that there is 
also a reaction between two different species Xj + X& — > X p , where j, k 7^ 1, the three 
dimensional set of equations for P(N\, Nj, Nk) should be included in order to keep track of 
the correlations between these species. The reaction network can be described by a graph in 
which the nodes represent the reactive species, and any two species that react are connected 
by an edge. In general, each set of multi-plane equations is associated with a maximal 
fully-connected subgraph, namely with a maximal set of nodes which are all connected to 
each other. Chemical networks tend to be sparse, namely most pairs of species do not react. 
Therefore, the multi-plane equations mostly consist of sets that invlove pairs of species and 
only few sets with three species or more. 

As an example, consider the case in which CO molecules are added to the H, O, OH 
system considered above . This gives rise to the following sequence of hydrogen addition 
reactions: H + CO -> HCO (X x + X 7 -> X 8 ), H + HCO -> H 2 CO {X x + X 8 -> X 9 ), H + 
H 2 CO -> H3CO (X x +X g -> X 10 ) and H + H 3 CO -> CH 3 OH {X x + X 10 -> X u ). Two other 
reactions that involve oxygen atoms also take place: O + CO — > C0 2 (X 2 + X 7 — >• X 12 ) and 
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O + HCO -> C0 2 + H (X 2 + X 8 -> X12 + Xi). This reaction network is described by the 
graph shown in Fig. 2. To account for the correlations between the populations of all pairs 
of species that react with each other, the multi-plane approach yields five sets of equations 
that account for P(N U N 3 ), P(N U N 9 ), P(N u N l0 ), P(Aq, 7V 2 , 7V 7 ) and P(N U N 2 ,N 8 ). The 
first three sets reside on planes while the last two sets are three dimensional sets that are 
needed in order to account for the correlations between the species that form CO2. Since the 
production of CO2 is low compared to other products, one may consider an approximation 
in which these correlations are neglected, reducing the equations into only planar sets. In 
this case the set of equations will include: P(N 1 ,N 2 ), P(N U N 3 ), P(N U N 7 ), P(Aq,7V 8 ), 
P(N U N 9 ) and P(iV\,iV 10 ). 

In the simulations presented here the activation energies for diffusion and desorption of 
CO, HCO, H 2 CO and H 3 CO were taken as E {7) = 30, #i(7) = 36, E (8) = 30, £1(8) = 38, 
E (9) = 33, £i(9) = 39 and E {10) = 35, Pi(10) = 41 meV, respectively. The flux of 
CO molecules was taken as F 7 = O.OOliq. In Fig. 3 we present the production rates per 
grain of H 2 0, C0 2 and CH 3 OH vs. grain size. An excellent agreement is found between the 
multi-plane approach (circles) and the complete master equation set (solid line). For small 
grains the approximated multi-plane set (+) shows significant deviations in the production 
rate of CO2 but still provides good results for the other species. The rate equation results 
are also shown (dashed line). For this reaction network the number of equations required in 
the complete master equation is between tens of thousands to several millions, depending on 
the grain size. Numerical integration of the complete master equation is impractical. The 
steady-state results presented here for the complete master equation set were obtained using 
Newton's method. Unlike the complete set, the multi-plane set consists of up to about a 
thousand equations, and can be simulated efficiently by numerical integration both under 
steady state and time dependent conditions using standard steppers such as Runge-Kutta. 

In summary, we have introduced the multi-plane method for the simulation of the master 
equation of grain-surface chemistry. For complex reaction networks, this method achieves a 
great reduction in the number of equations. It thus enables the incorporation of the master 
equation in models of interstellar chemistry, providing accurate results for the production 
rates of molecules on interstellar dust grains. We expect that the multi-plane method will 
also be useful in other problems that exhibit a related mathematical structure, such as 
the modeling of stratospheric cloud chemistry and the analysis of genetic networks in cells 



7 



2ll|22|. 

This work was supported by the Adler Foundation for Space Research of the Israel Science 
Foundation. 



[1] T.W. Hartquist and D.A. Williams, The chemically controlled cosmos (Cambridge University 

Press, Cambridge, UK, 1995). 
[2] R.J. Gould and E.E. Salpeter, Astropbys. J. 138, 393 (1963). 

[3] D. Hollenbach and E.E. Salpeter, J. Chem. Phys. 53, 79 (1970); Astrophys. J. 163, 155 (1971); 

D. Hollenbach, M.W. Werner and E.E. Salpeter, Astrophys. J. 163, 165 (1971). 
[4] J.B. Pickles and D.A. Williams, Astrophys. and Space Science 52, 433 (1977). 
[5] L.B. d'Hendecourt, L.J. Allamandola and J.M. Greenberg, Astron. Astrophys. 152, 130 

(1985). 

[6] PD. Brown and S.B. Charnley, Mon. Not. R. Astron. Soc. 244, 432 (1990). 

[7] T.I. Hasegawa and E. Herbst and CM. Leung, Astrophys. J. Supplement 82, 167 (1992). 

[8] P. Caselli, T.I. Hasegawa and E. Herbst, Astrophys. J. 408, 548 (1993). 

[9] A.G.G.M. Tielens and W. Hagen, Astron. Astrophys. 114, 245 (1982). 
[10] S.B. Charnley, A.G.G.M. Tielens and S.D. Rodgers, Astrophys. J. 482, L203 (1997). 
[11] P. Caselli, T.I. Hasegawa and E. Herbst, Astrophys. J. 495, 309 (1998). 
[12] O.M. Shalabiea, P. Caselli and E. Herbst, Astrophys. J. 502, 652 (1998). 
[13] T. Stantcheva, P. Caselli and E. Herbst, Astron. Astrophys. 375, 673 (2001). 
[14] O. Biham, I. Furman, V. Pirronello and G. Vidali, Astrophys. J. 553, 595 (2001). 
[15] N.J.B. Green, T. Toniazzo, M.J. Pilling, D.P. Ruffle, N. Bell and T.W. Hartquist, Astron. 

Astrophys. 375, 1111 (2001). 
[16] T. Stantcheva, V.I. Shematovich and E. Herbst, Astron. Astrophys. 391, 1069 (2002). 
[17] T. Stantcheva and E. Herbst, Mon. Not. R. Astron. Soc. 340, 983 (2003); Astron. Astrophys., 

in press (2004). 
[18] S.B. Charnley, Astrophys. J. 562, L99 (2001). 

[19] N. Katz, I. Furman, O. Biham, V. Pirronello and G. Vidali, Astrophys. J. 522, 305 (1999); 

H.B. Perets, O. Biham, V. Pirronello, J. Roser, S. Swords and G. Vidali, in preparation (2004). 
[20] J. Krug, Phys. Rev. E 67, 065102 (2003). 

8 



[21] H.H. McAdams and A. Arkin, Proc. Natl. Acad. Sci. USA 94, 814 (1997). 
[22] J. Paulsson and M. Ehrenberg, Phys. Rev. Lett. 84, 5447 (2000). 



9 



FIG. 1: The production rates of H2 (circles), O2 (squares), H2O (triangles) and the desorption rate 
of OH (x) per grain vs. grain size, given by the number of adsorption sites S and the diameter d. 
Excellent agreement is found between the multi-plane method (symbols) and the complete master 
equation (solid lines). The results of the rate equations are also shown (dashed lines). 

FIG. 2: A graph describing the reaction network that results from the adsorption of H and O 
atoms and CO molecules. The nodes represent the reactive species, while the edges connect pairs 
of species that react with each other. The reaction products are specified near the edges. 

FIG. 3: The production rates of H2O, CO2 and CH3OH vs. grain size as obtained from the 
multi-plane approach (circles) and the complete master equation (solid lines) are found to be in 
excellent agreement. The approximate multi-plane equations (+) deviate in the production rate of 
CO2, but are still accurate for the other species. The rate equation results are also shown (dashed 
lines) . 
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Fig. 3 



